

clear


%% data
data=readtable('../../../Data/data_main_sample.csv');
A=readtable('../../../Data/ddcg.csv'); % Acemoglu et al.'s (2019) replication data
CC=readtable('../../../Data/country_codes.csv');
% use 817 for Vietnam to merge
data.ccode(data.ccode==816)=817;

d1=data(data.year>=1875&data.year<=2000&isfinite(data.y)&isfinite(data.D),:);
d1(d1.year>=1914&d1.year<=1918,:)=[];
d1(d1.year>=1939&d1.year<=1945,:)=[];
d2=A(A.year>=1875&A.year<=2000&~strcmp(A.dem,'NA'),:);
d2(d2.year>=1914&d2.year<=1918,:)=[];
d2(d2.year>=1939&d2.year<=1945,:)=[];

ccodes=unique(d1.ccode);
years=unique(d1.year);
N=length(ccodes);
T=length(years);
y=zeros(N,T); Dbmr=y; D=y; M=y;
for n=1:N
    dn=d1(d1.ccode==ccodes(n),:);
    d2n=d2(strcmp(d2.wbcode,CC.wb{strcmp(CC.cown,num2str(ccodes(n)))}),:);
    for t=1:T
        if sum(years(t)==dn.year)>0
            if isfinite(dn.y(years(t)==dn.year))&&isfinite(dn.D(years(t)==dn.year))
                M(n,t)=1; % not missing
                y(n,t)=dn.y(years(t)==dn.year);
                Dbmr(n,t)=dn.D(years(t)==dn.year);
                if sum(years(t)==d2n.year)>0
                    D(n,t)=str2double(d2n.dem(years(t)==d2n.year));
                else
                    D(n,t)=Dbmr(n,t);
                end
            end
        end
    end
end
clear d1 dn n t A CC d2 d2n 

dd=readtable('../../../Data/distance_capitals.csv');
Z=10^(-6)*table2array(dd(:,2:end));


%% for Knitro
options=optimset('Display','off');


%% starting values

% "short" replication
load('tDGProbustDM_replication.mat','DGP')
pars0=DGP;

% % "long" replication
% load('tDGProbustDM_replication.mat','srng')
% rng(srng)
% feas=0; % ensure finite log_posterior starting value
% while feas==0
%     beta=[randn(N,1);randn(N,1)]/100;
%     s=rand(N,1);
%     ab=2*rand;
%     lambda=2*rand;
%     a=2*rand;
%     xi=2*rand;
%     b0A=randn/100;
%     b0D=randn/100;
%     pars0=[beta;s;ab;lambda;a;xi;b0A;b0D];
%     feas=isfinite(log_posterior_trueDGP(pars0,y,D,M,Z));
% end 


%% parameter estimates
[DGP,lp,exit]=knitromatlab(@(x)log_posterior_trueDGP(x,y,D,M,Z),pars0,[],[],...
  [],[],[-Inf(2*N,1);zeros(N+4,1);-Inf(2,1)],Inf(length(pars0),1),[],[],...
  options,'knitro.opt');


%% growth shocks variance matrix
load('../../model_prior/prior.mat')

s=DGP(2*N+(1:N));
a=DGP(3*N+3);
xi=DGP(3*N+4);

Sigma=diag(s)*(a*exp(-xi*Z))*diag(s);
Sigma=tril(Sigma)+tril(Sigma,-1)';  % enforce exact symmetry

clearvars -except a0 b0A b0D d_v d_vs f0 om_a om_b om_f om_th ...
    s_v s_vs th0 Sigma

save('priorWSigma_robustDM.mat')



